Optosurf simulation

Loading BokehJS ...

1 Gaussian function definition

Code
# Create a function to plot the equation
def plot_equation(mu, sigma, n, number_points, degrees, plot, title="Super-Gaussian", width = 700, height = 550):
    """
    Plot the Super-Gaussian equation using Bokeh

    Parameters
    ----------
    mu (float): Mean value for the equation
    sigma (float): Standard deviation value for the equation
    n (float): Order value for the equation
    number_points (int): number of points to calculate the function
    
    Returns
    -------
    plots (bokeh plot): Plot of the Super-Gaussian equation
    x(np): linspace for the gaussian plot
    y(np): gaussian values
    """
    # 1. Define linear degrees vector and calculate Super-Gaussian
    ticker = SingleIntervalTicker(interval=2.5, num_minor_ticks=10)
    xaxis = LinearAxis(ticker = ticker)
    x = np.linspace(degrees[0], degrees[1], number_points)
    y = np.exp(-abs(((x-mu)/sigma))**n)
    
    # 2. Plot 
    if plot:
        TOOLTIPS = [("index", "$index"),("(x,y)", "($x, $y)")]
        p = figure(title=title, x_axis_label='x', y_axis_label='y', tooltips = TOOLTIPS,
            width = width, height = height)
        p.line(x, y, line_width=4, alpha = 0.5, line_color = "#C5E064")
        p.add_layout(Grid(dimension=0, ticker=xaxis.ticker))
        p = plot_format(p, "Degrees", "Intensity", "bottom_left", "10pt", "10pt", "10pt")
        return p, x, y
    else:
        return x, y

p, x, y = plot_equation(0, 1, 3.5, 50000, [-15, 15], True)
show(p)

1.1 Window integration

The optosurf head has a 32 pixels line detectors. This is simulated by performing a window integration:

Code
# Create a function to do the window integration
def window_integration(number_windows, window_size, gap, x, y, p=None):
    """
    Performs a window integration

    Parameters
    ----------
    number_windows (int): Number of integration windows
    window_size (int): Number of data points in the window
    x(np): linspace for the gaussian plot
    y(np): gaussian values
    Returns
    -------
    p (bokeh plot): Plot of the integration
    integration_axis (np): window integration axis
    integration_points (np): Integrated points
    """
    integration_points = []
    integration_axis = []
    color_multiplier = len(bokeh.palettes.Viridis256)//number_windows
    count = 0
    
    for i in range(number_windows):
    # 1. Get data in every window and integrate
        a = i*window_size
        b = i*window_size + window_size
        
        x_temp = x[a:b-gap:1]
        y_temp = y[a:b-gap:1]
        integration = np.trapz(y_temp, x_temp, dx = x[1] - x[0])
        integration_points.append(integration)

        axis = x_temp[len(x_temp)//2]
        integration_axis.append(axis)

        # 2. Draw a rectangle of the integration window
        if p is not None:
            left_edge = x_temp[0]
            right_edge = x_temp[-1]
            p.rect(x=(left_edge + right_edge)/2, y=0.18, width=right_edge-left_edge, height=0.3, fill_alpha=0.001, fill_color='#C5E0B4', color='#C5E0B4')
            p.rect(x=(right_edge + x[b-1])/2, y=0.18, width=x[b-1]-right_edge, height=0.3, fill_alpha=0.005, fill_color='#F16C08', color = '#F16C08')
            p.circle(x_temp[::15], y_temp[::15], size = 4, alpha = 1)
            count += 1
    if p is not None:
        p.circle(integration_axis, integration_points, size = 7, color = '#FAA0A0')
        p.line(integration_axis, integration_points, line_width = 4, color = '#FAA0A0', alpha = 0.8)
    return p, integration_axis, integration_points

p, int_axis, int_points = window_integration(32, 1562, 100, x, y, p)
show(p)

The optosurf head has a 32 pixels line detectors. This is simulated by performing a window integration:

2 Histogram

The next step is to make a histogram reconstruction to calculate the standard deviation

Code
# Create a function to do histogram reconstruction
def histogram_reconstruction(int_points, hist_bool):
    """
    Constructs a histrogram

    Parameters
    ----------
    int_points(np): Points calculated from the window integration
    Returns
    -------
    hist_plot (bokeh plot): Plot of the histogram
    std_dev(float): Histogram's standard deviation
    """
    
    # a. Histogram reconstruction
    normalized_y = np.multiply(int_points, 10000)
    hist_2d = np.array([])
    for i, int_point in enumerate(normalized_y):
        round_int_point = round(float(int_point))
        hist_2d = np.concatenate((hist_2d, np.array([float(i)]*round_int_point)))
    
    # b. Calculate standard deviation
    stddev = np.std(hist_2d)
    
    # c. Plot histogram
    if hist_bool:
        hist_plot = alt.Chart(pd.DataFrame({'x': hist_2d})).mark_bar(opacity=0.7, color='#2D908C').encode(
        alt.X('x:Q', bin=alt.Bin(maxbins=20)),
        y='count()', 
        tooltip=['x','count()'])
        title = f'Standard deviation: {stddev}'
        hist_plot = hist_plot.properties(title=title, background='#282B30')
        return hist_plot, np.std(stddev)

    else:
        return stddev

alt.data_transformers.disable_max_rows()
alt.renderers.set_embed_options(theme='dark')
hist_plot, std_dev = histogram_reconstruction(int_points, True)
hist_plot

3 Standard deviation matrix

Code
mu_np = np.linspace(-5, 5, 4)
std_np = np.linspace(1.0, 5, 7)
X, Y = np.meshgrid(mu_np, std_np)
std_grid = np.empty_like(X)
plots_gaussian = []
n = 3.5
number_points = 10000
degrees = [-15, 15]
number_windows = 32
window_size = number_points//number_windows
gap = 100
for i in range(len(mu_np)):
    for j in range(len(std_np)):
        # Generate x and y Gaussian data points 
        title = f"mu: {mu_np[i]:.3f}, std: {std_np[j]:.3f}"
        # print(title)
        plot, x, y = plot_equation(mu_np[i], std_np[j], n, number_points, degrees, True, title, 200, 200)
        plot.title.text_font_size = "10pt"
        plot, int_axis, int_points = window_integration(number_windows, window_size, gap, x, y, plot)
        plots_gaussian.append(plot)

grid_gaussian = gridplot(children = plots_gaussian, ncols = len(std_np), merge_tools=False)
show(grid_gaussian)
Code
number_points = 5000
number_windows = 32
window_size = number_points//number_windows
gap = 100
n = 3.5
degrees = [-15, 15]
           
mu_range = [-5, 5]
axis_range = [mu_range[0], mu_range[1], 0.5]
axis_range = np.arange(-5, 5.5, 0.5)
mu_step = 0.10
mu_points = (mu_range[1] - mu_range[0])//mu_step

std_range = [1.0, 1.1]
std_step = 0.1 
std_points = (std_range[1] - std_range[0])/std_step


mu_np = np.linspace(mu_range[0], mu_range[1], int(mu_points))
std_np = np.linspace(std_range[0], std_range[1], int(std_points))
X, Y = np.meshgrid(mu_np, std_np)
std_grid = np.empty_like(X)

# Iterates mu and standard deviation
for i in range(len(mu_np)):
    for j in range(len(std_np)):
        # Generate x and y Gaussian data points 
        title = f"mu: {mu_np[i]:.3f}, std: {std_np[j]:.3f}"
        # print(title)
        x, y = plot_equation(mu_np[i], std_np[j], n, number_points, degrees, False)
        p, int_axis, int_points = window_integration(number_windows, window_size, gap, x, y)
        std_grid[j, i] = histogram_reconstruction(int_points, False)
        
source = pd.DataFrame({'mu': X.ravel(),
                        'std': Y.ravel(),
                        'value': std_grid.ravel()})

intensity_plot = alt.Chart(source).mark_rect().encode(
    alt.X('mu:O', axis=alt.Axis(values=axis_range, format=".1f")),
    y='std:O',
    color='value:Q',
    tooltip='value').interactive()
intensity_plot = intensity_plot.properties(width = 600, height = 200, title = 'Standard deviation variation')
intensity_plot